clear
clc
set(0,'DefaultFigureWindowStyle','docked')
addpath(genpath('library.model.20220926'))

%% about
% fits the following parameters:
% b, kappa, eta, sigma_z, 
% x_bar, sigma_x, chi
% to target the following moments:
% E[U], std[U], std[N*w]/std[Y], std[output growth vol], 
% E[phi|expansion], std[mu3_1y], std[1 year nom yld]

% load baseline (will use baseline process for phi)

load('run1a_baseline_sr.mat','parms')

sigma_x_base = parms.sigma_x;
rho_x_base = parms.rho_x;
x_bar_base = parms.x_bar;

clear parms

%% parameters

% grid sizes
parms.NN = 50;
parms.Nx = 50;
parms.Nz = 5;
% parms.Nz = 7;

% define recession threshold (actual recession prob depends on grid)
parms.RECESSION_PROB_THRESH = 1/6;

% grid: N
N_lb = 0.85; 
N_ub = 0.999;
parms.grid_N = linspace(N_lb, N_ub, parms.NN)';

% grid: z
parms.rho_z = 0.9; % 0.9: recessions (booms) last 12 (60) months on average
parms.sigma_z = 0.02; 

[parms.grid_z,parms.StateTransitionProbs] = rouwenhorst(0, parms.rho_z, parms.sigma_z^2, parms.Nz);

V = eig(parms.StateTransitionProbs');
cdf_z = cumsum(V(:,1)/sum(V(:,1)));

[~,parms.IDX_RECESSION_THRESH] = min(abs(cdf_z - parms.RECESSION_PROB_THRESH));

fprintf('Nz=%g, Prob(Recession)=%g (target is %g).\n', ...
    parms.Nz, cdf_z(parms.IDX_RECESSION_THRESH), parms.RECESSION_PROB_THRESH);

% labor market
parms.b = 0.81; % Unemployment benefit is b
parms.beta = 0.9979; %  Rep. household's time preference parameter 
parms.iota = 1.24; % Matching function
parms.phi_match = exp(-0.1/12); % constraint on matching (slow recovery)
% parms.phi_match = NaN; % constraint on matching (slow recovery), NaN=off
parms.kappa = 0.76; % vacancy costs
parms.eta = 0.158; % Bargaining power of worker 
parms.s = 0.034*ones(parms.Nz,1); % Hall AER (2017) following Shimer's calibration uses 0.0345

% preferences
parms.chi = 1/3.2; % CES aggregator
parms.gamma = 2; % risk aversion (intertemporal)

PHI_MAX = 0.99;
PHI_MIN = 0.1;
PHI_BAR = 0.9;

phi_grid = linspace(PHI_MIN,PHI_MAX,parms.Nx)';

parms.x_bar = -log(1/PHI_BAR - 1);
% parms.rho_x = 0.5^(1/12);
parms.rho_x = 0.9^(1/12);
% parms.rho_x = 0.85^(1/12);
parms.grid_x = -log(1./phi_grid - 1);
parms.sigma_x = 0.13*ones(parms.NN,parms.Nz); % should have sigma>0

parms.OPTION_x_standardize_shock = true; % dx loads on z'-E[z']/(std(z'))
% parms.OPTION_x_standardize_shock = false; % dx loads on z'-E[z']

%% check parameters

CheckParams_Model_20220926(parms)

%% Fit model

parms.x_bar = x_bar_base;
parms.rho_x = rho_x_base;
parms.sigma_x(:) = sigma_x_base(1); % should have sigma>0

NumTries = 1;
opts = optimset('display','iter-detailed','UseParallel',true,'MaxFunEvals',300);

U_mean_target = 0.0609;
U_std_target = 0.0078;
wagebillratio_target = 0.87;
vol_gy_target = 0.0081;
nom_yld_1y_mean_target = 5.28;
nom_yld_1y_vol_target = 3.32;
ave_mu3_vol_target = 0.0304*0.8^3;

kappa_guess = 0.065;
b_guess = 0.95;
eta_guess = 0.35;
sigma_z_guess = 0.0063;
chi_guess = 0.15;
beta_guess = 0.9979;
beta_lb = 0.99;
beta_ub = 0.9985;
log_sigma_x_guess = log(0.14);

% [parms, gesol, FitVal, FLAG_SUCCESSFUL_FSOLVE] = ...
%     FitModel20220926_Nstate_ver3(parms, kappa_guess, b_guess, eta_guess, sigma_z_guess, chi_guess, log_sigma_x_guess, ...
%         U_mean_target, U_std_target, wagebillratio_target, vol_gy_target, nom_yld_1y_vol_target, ave_mu3_vol_target, NumTries, opts);

[parms, gesol, FitVal, FLAG_SUCCESSFUL_FSOLVE] = ...
    FitModel20220926_Nstate_ver3(parms, kappa_guess, b_guess, eta_guess, sigma_z_guess, chi_guess, beta_guess, beta_lb, beta_ub, log_sigma_x_guess, ...
        U_mean_target, U_std_target, wagebillratio_target, vol_gy_target, nom_yld_1y_mean_target, nom_yld_1y_vol_target, ave_mu3_vol_target, NumTries, opts);

%% Verify solution
tic
Verify_Solution_Model_20220926(parms,gesol)
toc

%% compute model moments

[gesol,Phi_Nxz,prob_Nxz,Phi_NxYz,grid_Y,prob_NxYz,stats] = ...
    ComputeModelMoments_Model20220926_fast(parms,gesol);

%% save

save([mfilename,'.mat'],'parms','gesol','stats','FitVal', 'FLAG_SUCCESSFUL_FSOLVE')